import numpy as np
import matplotlib.pyplot as plt
import os

xscale = np.concatenate((np.array([3e10, 5e10, 5e3, 50, 200, 50, 150, 3000, 2e10, 3e10,5000, 100,
                   250, 100, 200, 3000, 37.0]), np.array([3, 200]), np.ones(1951)))
uscale = np.concatenate((np.array([3000/60, 3000/60, 3000/60, 3000/60, 2000, 2000, 50]), np.array([3000/60]), np.ones(1)))


class UtilsHelper:
    def load_ss(self, res_dir='data'):
        xss = np.load(os.path.join(res_dir, 'xss.npy'))
        uss = np.load(os.path.join(res_dir, 'uss.npy'))
        return xss, uss

    def modify_ss(self, xss, uss):
        """ If the modification is not required, simply return xss and uss as they are """
        c_inss = xss[14]  # Inlet concentration of mAb (mg/L)
        F_inss = uss[3]  # Inlet flow rate of buffer tank (L/min)
        # Steady state values of buffer tank are determined by the upstream
        F_outss = F_inss  # Outlet flow rate of buffer tank (L/min)
        xss_b = np.array([15, c_inss])  # Liquid height (dm), concentration of mAb (mg/L)
        uss_b = np.array([F_outss])
        # Steady state values for the integrated model
        xss_modify = np.concatenate((xss, xss_b, np.zeros(1951)))  # Downstream SS is 0, since it does not have SS
        uss_modify = np.concatenate((uss, uss_b, np.zeros(1)))  # Downstream SS is 0, since it does not have SS
        return xss_modify, uss_modify

    def prepare_ss(self, res_dir='data'):
        # Load steady state values of upstream
        xss, uss = np.array([0.483, 350.970]), np.array([299.709])  # With the computed economic zone
        # xss, uss = np.array((0.465, 352.000)), np.array((299.413))  # Without the computed economic zone
        # # Modify the steady state values. Project specific. Not necessarily useful for other projects
        # xss_modify, uss_modify = self.modify_ss(xss, uss)
        return xss, uss

    def prepare_init(self, xss, uss):
        """ Modify this function based on the demand """
        # Initialize lists for data storage
        t = []
        Xm = []  # MPC
        Um = []
        Xe = []  # EMPC
        Ue = []
        Xs = []  # Steady state
        Us = []
        Xi = []  # Plant integration
        Xie = []
        # Initialize states and inputs
        x0 = xss/xscale*0.95
        u0 = uss/uscale
        Xm += [x0]  # Scaled
        Xe += [x0]  # Scaled
        Xs += [xss] # Unscaled
        Xi += [x0]
        Xie += [x0]
        return Xm, Um, Xe, Ue, Xs, Us, Xi, Xie, u0, t

    def prepare_bounds(self, xss, uss):
        xlb = xss * 0.8
        xlb[16] = 33.0  # Temperature
        xlb[17] = 0  # Liquid level cannot be below the 0
        xlb[18] = xlb[14]  # [mAb] in buffer tank should be similar to that in the separator or reactor
        xlb[19:] = 0  # Prepare the lower bound of state for capture column

        xub = 1.2 * xss
        xub[7] = xscale[7]  # Volume of reactor
        xub[15] = xscale[15]  # Volume of separator
        xub[16] = xscale[16]  # Temperature
        xub[17] = 50  # For now, assume a big buffer tank
        xub[18] = xub[14]  # [mAb] in buffer tank should be similar to that in the separator or reactor
        for i in range(19,1969):  # (1950)
            if (i-19)%13 == 11:
                xub[i] = 36450  # qmax1
            elif (i-19)%13 == 12:
                xub[i] = 77850  # qmax2
            else:
                xub[i] = xub[14]  # [mAb] in capture column is similar to that in the separator
        xub[1969] = 60000  # Accumulated mAb

        ulb = 0.8 * uss
        ulb[-1] = int(0)
        uub = uscale
        uub[-1] = int(1)

        np.save('xlb', xlb)
        np.save('xub', xub)
        np.save('ulb', ulb)
        np.save('uub', uub)

        return xlb, xub, ulb, uub

    def plot_results(self, Xm, Um, Xe, Ue, Xs, Us, t, num_sim, dt_spl, cl):
        if cl:
            # convert results to numpy arrays
            t += [num_sim * dt_spl]
            Xm = np.array(Xm) * xscale
            Um = np.array(Um) * uscale
            Xe = np.array(Xe) * xscale
            Ue = np.array(Ue) * uscale
            Xs = np.array(Xs)
            Us = np.array(Us)
            self.compare_controllers(t, Xm, Um, Xe, Ue, Xs, Us)
            return Xm, Um, Xe, Ue, Xs, Us, t

        else:
            # convert results to numpy arrays
            t += [num_sim * dt_spl]
            Xm = np.array(Xm) * xscale
            Um = np.array(Um) * uscale
            Xs = np.array(Xs)
            Us = np.array(Us)
            X, U = self.process_results(Xm, Um)  # Convert results to numpy arrays
            # utils.save_results(t, X, U)  # Save results
            self.visualize_results_all(t, X, U)  # Visualize results
            return Xm, Um, Xe, Ue, Xs, Us, t

    def compare_controllers(self, t, Xm, Um, Xe, Ue, Xs, Us):
        t2 = t[:-1]  # Time list for inputs
        # create figure (fig), and array of axes (ax)
        plt.close("all")

        # Upstream
        plt.figure(1)
        plt.subplot(211)
        plt.plot(t, Xm[:, 0])
        plt.plot(t, Xe[:, 0])
        plt.plot(t, Xs[:, 0], linestyle="--")
        plt.ylabel('$X_{v, reactor}$ (cell/L)')
        plt.legend(['MPC', 'EMPC', 'SS'])
        plt.tight_layout()

        plt.subplot(212)
        plt.plot(t, Xm[:, 8])
        plt.plot(t, Xe[:, 8])
        plt.plot(t, Xs[:, 8], linestyle="--")
        plt.xlabel('Time (min)')
        plt.ylabel('$X_{v, separator}$ (cell/L)')
        plt.tight_layout()
        # plt.savefig("results/integratedmodel/upstream_state1.png")
        plt.show()

        plt.figure(2)
        plt.subplot(211)
        plt.plot(t, Xm[:, 1])
        plt.plot(t, Xe[:, 1])
        plt.plot(t, Xs[:, 1], linestyle="--")
        plt.ylabel('$X_{t, reactor}$ (cell/L)')
        plt.legend(['MPC', 'EMPC', 'SS'])
        plt.tight_layout()

        plt.subplot(212)
        plt.plot(t, Xm[:, 9])
        plt.plot(t, Xe[:, 9])
        plt.plot(t, Xs[:, 9], linestyle="--")
        plt.xlabel('Time (min)')
        plt.ylabel('$X_{t, separator}$ (cell/L)')
        plt.tight_layout()
        # plt.savefig("results/integratedmodel/upstream_state2.png")
        plt.show()

        plt.figure(3)
        plt.subplot(211)
        plt.plot(t, Xm[:, 2])
        plt.plot(t, Xe[:, 2])
        plt.plot(t, Xs[:, 2], linestyle="--")
        plt.ylabel('$[GLC]_{reactor}$ (mM)')
        plt.legend(['MPC', 'EMPC', 'SS'])
        plt.tight_layout()

        plt.subplot(212)
        plt.plot(t, Xm[:, 10])
        plt.plot(t, Xe[:, 10])
        plt.plot(t, Xs[:, 10], linestyle="--")
        plt.xlabel('Time (min)')
        plt.ylabel('$[GLC]_{separator}$ (mM)')
        plt.tight_layout()
        # plt.savefig("results/integratedmodel/upstream_state3.png")
        plt.show()

        plt.figure(4)
        plt.subplot(211)
        plt.plot(t, Xm[:, 3])
        plt.plot(t, Xe[:, 3])
        plt.plot(t, Xs[:, 3], linestyle="--")
        plt.ylabel('$[GLN]_{reactor}$ (mM)')
        plt.legend(['MPC', 'EMPC', 'SS'])
        plt.tight_layout()

        plt.subplot(212)
        plt.plot(t, Xm[:, 11])
        plt.plot(t, Xe[:, 11])
        plt.plot(t, Xs[:, 11], linestyle="--")
        plt.xlabel('Time (min)')
        plt.ylabel('$[GLN]_{separator}$ (mM)')
        plt.tight_layout()
        # plt.savefig("results/integratedmodel/upstream_state4.png")
        plt.show()

        plt.figure(5)
        plt.subplot(211)
        plt.plot(t, Xm[:, 4])
        plt.plot(t, Xe[:, 4])
        plt.plot(t, Xs[:, 4], linestyle="--")
        plt.ylabel('$[LAC]_{reactor}$ (mM)')
        plt.legend(['MPC', 'EMPC', 'SS'])
        plt.tight_layout()

        plt.subplot(212)
        plt.plot(t, Xm[:, 12])
        plt.plot(t, Xe[:, 12])
        plt.plot(t, Xs[:, 12], linestyle="--")
        plt.xlabel('Time (min)')
        plt.ylabel('$[LAC]_{separator}$ (mM)')
        plt.tight_layout()
        # plt.savefig("results/integratedmodel/upstream_state5.png")
        plt.show()

        plt.figure(6)
        plt.subplot(211)
        plt.plot(t, Xm[:, 5])
        plt.plot(t, Xe[:, 5])
        plt.plot(t, Xs[:, 5], linestyle="--")
        plt.xlabel('')
        plt.ylabel('$[AMM]_{reactor}$ (mM)')
        plt.legend(['MPC', 'EMPC', 'SS'])
        plt.tight_layout()

        plt.subplot(212)
        plt.plot(t, Xm[:, 13])
        plt.plot(t, Xe[:, 13])
        plt.plot(t, Xs[:, 13], linestyle="--")
        plt.xlabel('Time (min)')
        plt.ylabel('$[AMM]_{separator}$ (mM)')
        plt.tight_layout()
        # plt.savefig("results/integratedmodel/upstream_state6.png")
        plt.show()

        plt.figure(7)
        plt.subplot(211)
        plt.plot(t, Xm[:, 6])
        plt.plot(t, Xe[:, 6])
        plt.plot(t, Xs[:, 6], linestyle="--")
        plt.xlabel('')
        plt.ylabel('$[mAb]_{reactor}$ (mg/L)')
        plt.legend(['MPC', 'EMPC', 'SS'])
        plt.tight_layout()

        plt.subplot(212)
        plt.plot(t, Xm[:, 14])
        plt.plot(t, Xe[:, 14])
        plt.plot(t, Xs[:, 14], linestyle="--")
        plt.xlabel('Time (min)')
        plt.ylabel('$[mAb]_{separator}$ (mg/L)')
        plt.tight_layout()
        # plt.savefig("results/integratedmodel/upstream_state7.png")
        plt.show()

        plt.figure(8)
        plt.subplot(211)
        plt.plot(t, Xm[:, 7])
        plt.plot(t, Xe[:, 7])
        plt.plot(t, Xs[:, 7], linestyle="--")
        plt.xlabel('')
        plt.ylabel('$V_{reactor}$ (L)')
        plt.legend(['MPC', 'EMPC', 'SS'])
        plt.tight_layout()

        plt.subplot(212)
        plt.plot(t, Xm[:, 15])
        plt.plot(t, Xe[:, 15])
        plt.plot(t, Xs[:, 15], linestyle="--")
        plt.xlabel('Time (min)')
        plt.ylabel('$V_{separator}$ (L)')
        plt.tight_layout()
        # plt.savefig("results/integratedmodel/upstream_state8.png")
        plt.show()

        plt.figure(9)
        plt.subplot(211)
        plt.plot(t, Xm[:, 16])
        plt.plot(t, Xe[:, 16])
        plt.plot(t, Xs[:, 16], linestyle="--")
        plt.ylabel('T ($^\circ$C)')
        plt.legend(['MPC', 'EMPC', 'SS'])
        plt.tight_layout()

        plt.subplot(212)
        plt.step(t2, Um[:, 6], where='post')
        plt.step(t2, Ue[:, 6], where='post')
        plt.plot(t2, Us[:, 6], linestyle="--")
        plt.xlabel('Time (min)')
        plt.ylabel('$T_c$ ($^\circ$C)')
        plt.tight_layout()
        # plt.savefig("results/integratedmodel/upstream_state9.png")
        plt.show()

        plt.figure(10)
        plt.subplot(211)
        plt.step(t2, Um[:, 0], where='post')
        plt.step(t2, Ue[:, 0], where='post')
        plt.plot(t2, Us[:, 0], linestyle="--")
        plt.ylabel('$F_{in}$ (L/min)')
        plt.legend(['MPC', 'EMPC', 'SS'])
        plt.tight_layout()

        plt.subplot(212)
        plt.step(t2, Um[:, 1], where='post')
        plt.step(t2, Ue[:, 1], where='post')
        plt.plot(t2, Us[:, 1], linestyle="--")
        plt.xlabel('Time (min)')
        plt.ylabel('$F_1$ (L/min)')
        plt.tight_layout()
        # plt.savefig("results/integratedmodel/upstream_input1.png")
        plt.show()

        plt.figure(11)
        plt.subplot(211)
        plt.step(t2, Um[:, 2], where='post')
        plt.step(t2, Ue[:, 2], where='post')
        plt.plot(t2, Us[:, 2], linestyle="--")
        plt.ylabel('$F_r$ (L/min)')
        plt.legend(['MPC', 'EMPC', 'SS'])
        plt.tight_layout()

        plt.subplot(212)
        plt.step(t2, Um[:, 3], where='post')
        plt.step(t2, Ue[:, 3], where='post')
        plt.plot(t2, Us[:, 3], linestyle="--")
        plt.xlabel('Time (min)')
        plt.ylabel('$F_2$ (L/min)')
        plt.tight_layout()
        # plt.savefig("results/integratedmodel/upstream_input2.png")
        plt.show()

        plt.figure(12)
        plt.subplot(211)
        plt.step(t2, Um[:, 4], where='post')
        plt.step(t2, Ue[:, 4], where='post')
        plt.plot(t2, Us[:, 4], linestyle="--")
        plt.ylabel('$[GLC]_{in}$ (mM)')
        plt.legend(['MPC', 'EMPC', 'SS'])
        plt.tight_layout()

        plt.subplot(212)
        plt.step(t2, Um[:, 5], where='post')
        plt.step(t2, Ue[:, 5], where='post')
        plt.plot(t2, Us[:, 5], linestyle="--")
        plt.xlabel('Time (min)')
        plt.ylabel('$[GLN]_{in}$ (mM)')
        plt.tight_layout()
        # plt.savefig("results/integratedmodel/upstream_input3.png")
        plt.show()

        # Buffer tank
        plt.figure(13)
        plt.subplot(211)
        plt.step(t2, Um[:, 3], where='post')
        plt.step(t2, Ue[:, 3], where='post')
        plt.plot(t2, Us[:, 3], linestyle="--")
        plt.ylabel('$F_{in, buffer tank}$ (L/min)')
        plt.legend(['MPC', 'EMPC', 'SS'])
        plt.tight_layout()

        plt.subplot(212)
        plt.step(t2, Um[:, 7], where='post')
        plt.step(t2, Ue[:, 7], where='post')
        plt.plot(t2, Us[:, 7], linestyle='--')
        plt.xlabel('Time (min)')
        plt.ylabel('$F_{out, buffer tank}$ (L/min)')
        plt.legend(['PID-MPC', 'PID-EMPC', 'SS'])
        plt.tight_layout()
        # plt.savefig("results/integratedmodel/buffer_input1.png")
        plt.show()

        plt.figure(14)
        plt.subplot(211)
        plt.plot(t, Xm[:, 17])
        plt.plot(t, Xe[:, 17])
        plt.plot(t, Xs[:, 17], linestyle="--")
        plt.ylabel('$h_{buffertank}$ (dm)')
        plt.legend(['MPC', 'EMPC', 'SS'])
        plt.tight_layout()

        plt.subplot(212)
        plt.plot(t, Xm[:, 18])
        plt.plot(t, Xe[:, 18])
        plt.plot(t, Xs[:, 18], linestyle="--")
        plt.xlabel('Time (min)')
        plt.ylabel('$[mAb]_{buffertank}$ (mg/L)')
        plt.tight_layout()
        # plt.savefig("results/integratedmodel/buffer_state1.png")
        plt.show()

        # Downstream
        plt.figure(11)
        plt.subplot(231)
        plt.plot(t, Xm[:, 19+13 * 1 + 1])
        plt.plot(t, Xm[:, 19+13 * 75 + 1])
        plt.plot(t, Xm[:, 19+13 * 149 + 1])
        plt.plot(t, Xe[:, 19+13 * 1 + 1], linestyle='--')
        plt.plot(t, Xe[:, 19+13 * 75 + 1], linestyle='--')
        plt.plot(t, Xe[:, 19+13 * 149 + 1], linestyle='--')
        # plt.ylim((0.88, 0.91))
        plt.ylabel("$c_p$ (mg/L) r=1")
        plt.xlabel("Time (min)")
        plt.subplot(232)
        plt.plot(t, Xm[:, 19+13 * 1 + 5])
        plt.plot(t, Xm[:, 19+13 * 75 + 5])
        plt.plot(t, Xm[:, 19+13 * 149 + 5])
        plt.plot(t, Xe[:, 19+13 * 1 + 5], linestyle='--')
        plt.plot(t, Xe[:, 19+13 * 75 + 5], linestyle='--')
        plt.plot(t, Xe[:, 19+13 * 149 + 5], linestyle='--')
        # plt.ylim((0.88, 0.91))
        plt.ylabel("$c_p$ (mg/L) r=5")
        plt.xlabel("Time (min)")
        plt.subplot(233)
        plt.plot(t, Xm[:, 19+13 * 1 + 10])
        plt.plot(t, Xm[:, 19+13 * 75 + 10])
        plt.plot(t, Xm[:, 19+13 * 149 + 10])
        plt.plot(t, Xe[:, 19+13 * 1 + 10], linestyle='--')
        plt.plot(t, Xe[:, 19+13 * 75 + 10], linestyle='--')
        plt.plot(t, Xe[:, 19+13 * 149 + 10], linestyle='--')
        # plt.ylim((0.88, 0.91))
        plt.ylabel("$c_p$ (mg/L) r=10")
        plt.xlabel("Time (min)")
        plt.subplot(234)
        plt.plot(t, Xm[:, 19+13 * 1 + 0])
        plt.plot(t, Xm[:, 19+13 * 75 + 0])
        plt.plot(t, Xm[:, 19+13 * 149 + 0])
        plt.plot(t, Xe[:, 19+13 * 1 + 0], linestyle='--')
        plt.plot(t, Xe[:, 19+13 * 75 + 0], linestyle='--')
        plt.plot(t, Xe[:, 19+13 * 149 + 0], linestyle='--')
        # plt.ylim((0.88, 0.91))
        plt.ylabel("$[mAb]_{column}$ (mg/L)")
        plt.xlabel("Time (min)")
        plt.subplot(235)
        plt.plot(t, Xm[:, 19+13 * 1 + 11])
        plt.plot(t, Xm[:, 19+13 * 75 + 11])
        plt.plot(t, Xm[:, 19+13 * 149 + 11])
        plt.plot(t, Xe[:, 19+13 * 1 + 11], linestyle='--')
        plt.plot(t, Xe[:, 19+13 * 75 + 11], linestyle='--')
        plt.plot(t, Xe[:, 19+13 * 149 + 11], linestyle='--')
        # plt.ylim((32, 35))
        plt.ylabel("$q_1$ (mg/L)")
        plt.xlabel("Time (min)")
        plt.subplot(236)
        plt.plot(t, Xm[:, 19+13 * 1 + 12])
        plt.plot(t, Xm[:, 19+13 * 75 + 12])
        plt.plot(t, Xm[:, 19+13 * 149 + 12])
        plt.plot(t, Xe[:, 19+13 * 1 + 12], linestyle='--')
        plt.plot(t, Xe[:, 19+13 * 75 + 12], linestyle='--')
        plt.plot(t, Xe[:, 19+13 * 149 + 12], linestyle='--')
        # plt.ylim((71, 74))
        plt.ylabel("$q_2$ (mg/L)")
        plt.xlabel("Time (min)")
        plt.legend(['MPC level 1', 'level 75', 'level 149', 'EMPC level 1', 'level 75', 'level 149'])
        plt.tight_layout()
        # plt.savefig("results/integratedmodel/downstream_state1.png")
        plt.show()

    def visualize_results_all(self, t, X, U):
        t2 = t[:-1]  # Time list for inputs
        plt.close("all")

        # States
        plt.figure(1)
        plt.subplot(211)
        plt.plot(t, X[:, 0])
        plt.ylabel('$X_{v, reactor}$ (cell/L)')
        plt.subplot(212)
        plt.plot(t, X[:, 8])
        plt.xlabel('Time (min)')
        plt.ylabel('$X_{v, separator}$ (cell/L)')
        # plt.savefig("results/integratedmodel_ol/upstream_state1.png")
        plt.show()

        plt.figure(2)
        plt.subplot(211)
        plt.plot(t, X[:, 1])
        plt.ylabel('$X_{t, reactor}$ (cell/L)')
        plt.subplot(212)
        plt.plot(t, X[:, 9])
        plt.xlabel('Time (min)')
        plt.ylabel('$X_{t, separator}$ (cell/L)')
        # plt.savefig("results/integratedmodel_ol/upstream_state2.png")
        plt.show()

        plt.figure(3)
        plt.subplot(211)
        plt.plot(t, X[:, 2])
        plt.ylabel('$[GLC]_{reactor}}$ (mM)')
        plt.subplot(212)
        plt.plot(t, X[:, 10])
        plt.xlabel('Time (min)')
        plt.ylabel('$[GLC]_{separator}}$ (mM)')
        # plt.savefig("results/integratedmodel_ol/upstream_state3.png")
        plt.show()

        plt.figure(4)
        plt.subplot(211)
        plt.plot(t, X[:, 3])
        plt.ylabel('$[GLN]_{reactor}$ (mM)')
        plt.subplot(212)
        plt.plot(t, X[:, 11])
        plt.xlabel('Time (min)')
        plt.ylabel('$[GLN]_{separator}$ (mM)')
        # plt.savefig("results/integratedmodel_ol/upstream_state4.png")
        plt.show()

        plt.figure(5)
        plt.subplot(211)
        plt.plot(t, X[:, 4])
        plt.ylabel('$[LAC]_{reactor}$ (mM)')
        plt.subplot(212)
        plt.plot(t, X[:, 12])
        plt.xlabel('Time (min)')
        plt.ylabel('$[LAC]_{separator}$ (mM)')
        # plt.savefig("results/integratedmodel_ol/upstream_state5.png")
        plt.show()

        plt.figure(6)
        plt.subplot(211)
        plt.plot(t, X[:, 5])
        plt.ylabel('$[AMM]_{reactor}$ (mM)')
        plt.subplot(212)
        plt.plot(t, X[:, 13])
        plt.xlabel('Time (min)')
        plt.ylabel('$[AMM]_{separator}$ (mM)')
        # plt.savefig("results/integratedmodel_ol/upstream_state6.png")
        plt.show()

        plt.figure(7)
        plt.subplot(211)
        plt.plot(t, X[:, 6])
        plt.ylabel('$[mAb]_{reactor}$ (mg/L)')
        plt.subplot(212)
        plt.plot(t, X[:, 14])
        plt.xlabel('Time (min)')
        plt.ylabel('$[mAb]_{separator}$ (mg/L)')
        # plt.savefig("results/integratedmodel_ol/upstream_state7.png")
        plt.show()

        plt.figure(8)
        plt.subplot(211)
        plt.plot(t, X[:, 7])
        plt.ylabel('$V_{reactor}$ (L)')
        plt.subplot(212)
        plt.plot(t, X[:, 15])
        plt.xlabel('Time (min)')
        plt.ylabel('$V_{separator}$ (L)')
        # plt.savefig("results/integratedmodel_ol/upstream_state8.png")
        plt.show()

        plt.figure(9)
        plt.subplot(211)
        plt.plot(t, X[:, 16])
        plt.ylabel('T ($^\circ$C)')
        plt.subplot(212)
        plt.step(t2, U[:, 6], where='post')
        plt.xlabel('Time (min)')
        plt.ylabel('$T_c$ ($^\circ$C)')
        # plt.savefig("results/integratedmodel_ol/upstream_state9.png")
        plt.show()

        plt.figure(10)
        plt.subplot(211)
        plt.plot(t, X[:, 17])
        plt.ylabel('$h_{buffer tank}$ (dm)')
        plt.subplot(212)
        plt.plot(t, X[:, 18])
        plt.xlabel('Time (min)')
        plt.ylabel('$[mAb]_{buffer tank}$ (mg/L)')
        # plt.savefig("results/integratedmodel_ol/buffer_state1.png")
        plt.show()

        plt.figure(11)
        plt.subplot(231)
        plt.plot(t, X[:, 19+13 * 1 + 1])
        plt.plot(t, X[:, 19+13 * 75 + 1])
        plt.plot(t, X[:, 19+13 * 149 + 1])
        # plt.ylim((0.88, 0.91))
        plt.ylabel("$c_p$ (mg/L) r=1")
        plt.xlabel("t (min)")
        plt.subplot(232)
        plt.plot(t, X[:, 19+13 * 1 + 5])
        plt.plot(t, X[:, 19+13 * 75 + 5])
        plt.plot(t, X[:, 19+13 * 149 + 5])
        # plt.ylim((0.88, 0.91))
        plt.ylabel("$c_p$ (mg/L) r=5")
        plt.xlabel("t (min)")
        plt.subplot(233)
        plt.plot(t, X[:, 19+13 * 1 + 10])
        plt.plot(t, X[:, 19+13 * 75 + 10])
        plt.plot(t, X[:, 19+13 * 149 + 10])
        # plt.ylim((0.88, 0.91))
        plt.ylabel("$c_p$ (mg/L) r=10")
        plt.xlabel("t (min)")
        plt.subplot(234)
        plt.plot(t, X[:, 19+13 * 1 + 0])
        plt.plot(t, X[:, 19+13 * 75 + 0])
        plt.plot(t, X[:, 19+13 * 149 + 0])
        # plt.ylim((0.88, 0.91))
        plt.ylabel("$[mAb]_{column}$ (mg/L)")
        plt.xlabel("t (min)")
        plt.subplot(235)
        plt.plot(t, X[:, 19+13 * 1 + 11])
        plt.plot(t, X[:, 19+13 * 75 + 11])
        plt.plot(t, X[:, 19+13 * 149 + 11])
        # plt.ylim((32, 35))
        plt.ylabel("$q_1$ (mg/L)")
        plt.xlabel("t (min)")
        plt.subplot(236)
        plt.plot(t, X[:, 19+13 * 1 + 12])
        plt.plot(t, X[:, 19+13 * 75 + 12])
        plt.plot(t, X[:, 19+13 * 149 + 12])
        # plt.ylim((71, 74))
        plt.ylabel("$q_2$ (mg/L)")
        plt.xlabel("t (min)")
        plt.legend(['level 1', 'level 75', 'level 149'])
        plt.tight_layout()
        # plt.savefig("results/integratedmodel_ol/downstream_state1.png")
        plt.show()

        # Inputs
        plt.figure(12)
        plt.subplot(211)
        plt.step(t2, U[:, 0], where='post')
        plt.ylabel('$F_{in}$ (L/min)')
        plt.subplot(212)
        plt.step(t2, U[:, 1], where='post')
        plt.xlabel('Time (min)')
        plt.ylabel('$F_1$ (L/min)')
        # plt.savefig("results/integratedmodel_ol/upstream_input1.png")
        plt.show()

        plt.figure(13)
        plt.subplot(211)
        plt.step(t2, U[:, 2], where='post')
        plt.ylabel('$F_r$ (L/min)')
        plt.subplot(212)
        plt.step(t2, U[:, 3], where='post')
        plt.xlabel('Time (min)')
        plt.ylabel('$F_2$ (L/min)')
        # plt.savefig("results/integratedmodel_ol/upstream_input2.png")
        plt.show()

        plt.figure(14)
        plt.subplot(211)
        plt.step(t2, U[:, 4], where='post')
        plt.ylabel('$[GLC]_{in}$ (mM)')
        plt.subplot(212)
        plt.step(t2, U[:, 5], where='post')
        plt.xlabel('Time (min)')
        plt.ylabel('$[GLN]_{in}$ (mM)')
        # plt.savefig("results/integratedmodel_ol/upstream_input3.png")
        plt.show()

        plt.figure(16)
        plt.subplot(211)
        plt.step(t2, U[:, 3], where='post')
        plt.ylabel('$F_{in, buffer tank}$ (L/min)')
        plt.subplot(212)
        plt.step(t2, U[:, 7], where='post')
        plt.xlabel('Time (min)')
        plt.ylabel('$F_{out, buffer tank}$ (L/min)')
        # plt.savefig("results/integratedmodel_ol/buffer_input1.png")
        plt.show()

        t_special = [time_instant/60 for time_instant in t]
        plt.figure(17)
        plt.plot(t_special, X[:, 19+13 * 1 + 0])
        plt.plot(t_special, X[:, 19+13 * 75 + 0])
        plt.plot(t_special, X[:, 19+13 * 149 + 0])
        plt.ylabel("$[mAb]_{column}$ (mg/L)")
        plt.xlabel("t (hr)")
        plt.legend(['level 1', 'level 75', 'level 149'])
        plt.tight_layout()
        # plt.savefig("results/integratedmodel_ol/downstream_state1.png")
        plt.show()

    def process_results(self, X, U):
        X = np.array(X)
        U = np.array(U)
        return X, U

    def plot_learning_curve(self, scores, x, figure_file):
        running_avg = np.zeros(len(scores))
        for i in range(len(running_avg)):
            running_avg[i] = np.mean(scores[max(0, i - 100):(i + 1)])
        plt.plot(x, running_avg)
        plt.title('Running average of previous 100 scores')
        plt.savefig(figure_file)
        # plt.show()